A study of the Navier-Stokes-a model for two-dimensional turbulence 
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ABSTRACT. The Navier-Stokes-a sub-grid scale model of turbulence is a mollification of the Navier-Stokes 
equations in which the vorticity is advected and stretched by a smoothed velocity field. The smoothing is per- 
formed by filtering the velocity field over spatial scales of size smaller than a. This is achieved by convolution 
with a kernel associated with the Green's function of the Helmholtz operator scaled by a parameter a. The 
statistical properties of the smoothed velocity field are expected to match those of Navier-Stokes turbulence for 
scales larger than a, thus providing a more computable model for those scales. For wavenumbers k such that 
ka > 1, corresponding to spatial scales smaller than a, there are three candidate power laws for the energy 
spectrum, corresponding to three possible dynamical eddy turnover time scales in the model equations: one from 
the smoothed field, the second from the rough field and the third from a special combination of the two. Using 
two-dimensional turbulence as a test case, we measure the scaling of the spectra from high-resolution simulations 
of the Navier-Stokes-a model, in the limit as a— >co. We show that the energy spectrum of the smoothed velocity 
field scales as k~ 7 in the direct enstrophy cascade regime, consistent with dynamics dominated by the time scale 
of the rough velocity field. This result implies that the rough velocity field plays a crucial role in the development 
of the smoothed field even though the latter is the purported model for the Navier-Stokes turbulence. This effect 
must be taken into account when performing accurate simulations of the NS-a model. 
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Let v(x,t) denote the velocity field and p(x,t) the pressure of a homogeneous incom- 
pressible fluid with constant density p(x) = 1. The Navier-Stokes equations (NSE): 
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1. Introduction 




d t v — v x (V x v) = —Vp + uAv + /, 
V • v = 0, 

v(x,0) = v in (x), 
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govern the dynamics of such fluid flows, where the body force f(x), the viscosity v > 0, 
and the initial velocity field v m (x) are given. We supplement the system in (1) with 
periodic boundary conditions in a basic box [0, L] n , where n = 2, or 3. We assume here, 
f\ 0L ]nf( x ) dx = L L j n v m (x) dx = 0, which implies that L L , n v(x, t) dx = for all 
t > 0. It is computationally prohibitively expensive to perform reliable direct numerical 
simulation (DNS) of the NSE for high Reynolds number flows, due to the wide range of 
scales of motion that need to be resolved. The use of numerical models allows researchers 
to simulate turbulent flows of interest using smaller computational resources. In this 
paper we study a particular sub-grid scale turbulence model known as the Navier-Stokes- 
oi (NS-a) model: 

d t v — u x (V x v) = — Vp + vAv + / 
(2) V ■ u = V ■ v = 

v = u — a 2 Am. 

The system in (2) reduces to the Navier-Stokes system when a = 0, i.e. u = v. One 
can think of the parameter a as the length scale associated with the width of the filter 
which smooths v to obtain u. The filter is associated with the Green's function of the 
Helmholtz operator (J — a 2 A) (given in the third equation of (2)). When a > 0, notice 
that the nonlinear term u x (V x v) in (2) is milder than that in (1). 

The inviscid and unforced version of the a-model in (2) was derived in [19] based on the 
Hamilton variational principle subject to the incompressibility constraint div v = 0. It 
is worth mentioning that the inviscid model, the so called Euler-a model, coincides with 
the inviscid second-grade non-Newtonian fluid model (see, e.g., [8, 9, 12]). By adding the 
viscous term — uAv and the forcing / in an ad hoc fashion, the authors in [3, 4, 5] and [14] 
obtain the NS-a system (2) which they then named the viscous Camassa-Holm equations 
(VCHE), also known as the Lagrangian averaged Navier-Stokes-a model (LANS-a). The 
LANS-a model is different from the viscous second-grade non-Newtonian visco-elastic 
fluid in that the viscous term in the former is — uAv , distinguishable from — uAu in the 
latter. In references [3, 4, 5] it was found that the analytical steady state solutions for 
the 3d NS-Q! model compared well with averaged experimental data from turbulent flows 
in channels and pipes for wide range of Reynolds numbers. It was this fact which led the 
authors of [3, 4, 5] to suggest that the NS-a model be used as a closure model for the 
Reynolds averaged equations (RANS). It has since been discovered that there is in fact 
a whole family of 'a'- models which provide similar successful comparison with empirical 
data - among these are the Clark-a model [1, 10], the Leray-a model [7], the modified 
Leray-a model [20] and the simplified Bardina model [2, 24] (see also [27] for a family of 
similar models). 

We assume periodic boundary conditions with basic box [0, L] n , and denote the spatial 
Fourier coefficients of u{x) and v(x) by u(£) and v(£), respectively. Let (•, •) denote the 
L 2 -inner product and | • | denote the L 2 -norm over the basic box [0, L] n (we also use | • | 
for the modulus of a vector). In the two-dimensional (2d) case we have two conserved 



quantities for the (inviscid and unforced) NS-a equations, namely the energy, 

oo 

(3) E a = E <*(0 where > 

I5l=i 

(4) E a (t) = l(m-m) = \(\m\ 2 +*\\ 2 \m\ 2 ) 

and the enstrophy, 

oo 

(5) Q a = 5^0,(0 

151=1 

(6) where fi a (0 = ^(£)| 2 

(7) and the vorticity q(£) = (V x n)((). 
The energy of the smoothed field u is given by, 

oo 

(8) E u = J2 EU (0 

l*l=i 

(9) where £ u (£) = i|fi(£)| 2 , 

and is not conserved in the inviscid and unforced equations. 

For a wavenumber k, we define the component u k of a velocity field u by 

(10) u k := u k {x) = HO^ 

|£|=k 

and the component Uk>,k" with a range of wavenumbers [k', k") by 

(11) U k ',k" •= Uk',k"(x) = ^ Mfe. 

fc'<fc<fc" 

Denote by E a (k) the energy spectrum associated to 

(12) t% = ^(\u k , 2k \ 2 + a 2 \Vu k , 2k \ 2 ), 

which is the average energy per unit mass of eddies of linear size I E ( — -, — , then 

(13) t a k = / E a ( X ) d X - 

Jk 

Similarly, denote by E u (k) the energy spectrum associated to ej£ = - (|^fc,2fc| 2 ), that is, 

p2k 

(14) 4 = / E U { X ) d X . 
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In [15], it can be shown that E u (k) has the following remarkable property in three- 
dimensional (3d) viscous NS-a case, 



(15) E u {k) 



k 5 / 3 , when to < 1 
k~P, when ka ^> 1 



where — 3 > 5/3, corresponds to a sharper roll-off than the NSE scaling of the energy 
power spectrum k~ 5 ^ 3 for ka ^> 1. Thus, E u {k) is offered as a model spectrum for the NSE 
for spatial scales of motion larger than a. From the point of view of numerical simulation, 
the faster (compared to NSE) decay of E u {k) of NS-a for ka > 1, indicates suppression 
of scales smaller than a and in principle implies reduced resolution requirements for 
simulating a given flow. The scaling in the ka > 1 regime, that is, the exponent j3, was 
determined in [14] using a time scale which depends solely on the smoothed velocity field. 
Later on, in [1, 7, 20], three possibilities for the exponent (3 were derived since there 
are two velocities in the NS-a model. The three possibilities for the exponent j3 depend 
on whether the turnover time scale of an eddy of size k~ x is determined by (A;|wfc|) -1 , 
(klvkl)' 1 , or (ky/ (uk, v k ) ) _1 (see e.g. [1, 20, 2, 15] and Section 2). Determining the 
actual scaling requires resolved numerical simulations. In the rest of the paper we will 
use the notation E(k) and E u (k) interchangeably for finite a simulations when there is 
no ambiguity in meaning. 

The 3d NS-a model was tested numerically in [6], for moderate Reynolds number in 
a simulation of size 256 3 , with periodic boundary conditions. It was observed that the 
large scale features of a turbulent flow were indeed captured and there was a rollover of 
the energy spectrum from k~ 5 ^ 3 for ka <C 1 to something steeper for ka ^> 1, although 
the scaling ranges were insufficient to enable extraction of the scaling exponent f3 of (15) 
unambiguously. Other numerical tests of the NS-a model were performed in [17], [18], 
and [25], with similar results. 

Our goal in this work is to measure the scaling of the energy spectrum of the NS-a in 
the regime ka ^> 1 in the forward enstrophy cascade regime for 2d turbulence. Assuming 
the validity of the semi-rigorous arguments introduced in [15], we can then infer from our 
computations the timescale which governs the turnover of an eddy of size 1/k. We choose 
to do this measurement for the 2d case with the expectation that discerning scaling ranges 
and exponents will be cheaper and more tenable in a 2d system than in a 3d system at 
the same grid-resolution. Nevertheless, we hope that the results in 2d will provide some 
insight into the turnover time of eddies of size 1/k for ka 3> 1 in the 3d case. 

In two-dimensional turbulence there are two inertial ranges, one for the inverse cascade 
of energy E and the other for the direct cascade of enstrophy with corresponding scalings 
of the energy spectrum as follows [22]: 



E(k) 



k 5//3 , where k -C kf, inverse energy cascade regime 



k 3 , where kd ^> k ^> kf, direct enstrophy cascade regime 
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where kf is the injection (forcing) wavenumber for both energy and enstrophy, and kd is 
the enstrophy dissipation wavenumber. The k~ 3 scaling is thought to be accurate upto 
a log- correct ion [23]. In a NS-a model implementation, if we choose a smaller than the 
injection scale such that kjot <C 1, we expect a modification of the scaling in the enstrophy 
cascade regime as follows 

{/c~ 5 / 3 , where k <C 
k~ 3 , where k f < k < 1/a, > 
fc~ 7 , where 1 / a <ti k <ti k d J 

where 7 > 3. For the 2d NS-a, we show, using semi-rigorous arguments as in [1, 13, 20, 
2, 15] that there are three possible scalings, 7 G {7, 19/3, 17/3}. We verify numerically 
which of these scalings actually emerges, and therefore deduce the eddy turnover time for 
the eddies of size smaller than the length scale a. 

As may be seen from Eqs. (16), for finite a, there are three distinct scaling ranges for 
2d NS-a. Since we are primarily interested in the value of 7 we would like to maximize 
the range of scales over which A; -7 scaling dominates. Therefore, we minimize the inverse 
cascade range of k~ 5 ^ 3 scaling by forcing in the lowest modes. We also minimize the k~ 3 
scaling range by increasing a all the way to a — > 00. The latter limit is well-defined and 
yields what we will call the NS-00 equations (see, e.g. [21]): 

d t v — u x V x v = — Vp + vAv + / 

(17) V ■ u = V • v = 

v = -L 2 Au, 

where the forcing term was rescaled appropriately to avoid trivial dynamics. Assuming 
that the scaling of the spectrum as a — > 00 is identical to the scaling in the range 
ka 3> 1, for finite (small) a and sufficiently long scaling ranges, we obtain the first 
resolved numerical calculation of the sub-a scales. The assumption is not unreasonable as 
nothing singular is expected to occur as a grows, other than the lengthening of the scaling 
range of interest. We remark, however, that one has to rescale the forcing appropriately 
in order to avoid trivial dynamics for large values of a, and decaying turbulence at the 
limit when a— >oc. 

The paper is organized as follows. In section 2 we summarize our analytical re- 
sults. We show that the 2d NS-a model exhibits an inverse cascade of the energy 
E a = I Jt ^ u(x)-v(x) dx and a forward cascade of the enstrophy Q a = ~ Jr 0L i 2 k(a ; )| 2 dx, 
where q(x) = V x v(x). Then, we show that in the forward enstrophy cascade regime, 
the energy of the smoothed velocity E u — | Jj L j 2 u(x) ■ u(x) dx should follow the Kraich- 

nan k~ 3 power law for ka <C 1. For k such that ka ^> 1 and kf k kd, we will show 
how the three possible power laws, k~ 17 ^ 3 , k~ 19 ^ 3 and k~ 7 , arise. 

In section 3 we present the details of our numerical scheme, data parameters and 
empirical justification for adopting a hypoviscousity sink term for the energy in the low 
modes. Our numerical calculations ranged in resolution from 256 2 to 4096 2 . Our data 
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show a convergence to k~ 7 power law in the range ka <C 1 as a — > oo. This scaling 
allows us to conclude that the eddy turnover time in the enstrophy cascade region of the 
smoothed 2d NS-a energy spectrum is determined by the time scale which depends on 
the square of the unsmoothed velocity field, that is the variable v. 

2. Statistical predictions for the NS-a model 

2.1. Average transfer and cascade of energy and enstrophy for the 2d NS-a 

model. In two-dimensional turbulence, energy and enstrophy transfer behave as follows: 
in the wavenumber range below kf, the energy and enstrophy go from higher modes to 
lower modes; in the wavenumber range above kf, the energy and enstrophy go from low 
to high wavenumbers. In a certain wavenumber regime above kf, there is much stronger 
direct transfer of enstrophy than energy, leading to what is called the direct enstrophy 
cascade. Similarly, in a certain range below kf, there is a more dominant transfer of 
energy, leading to what is called the inverse energy cascade. For the 2d NS-a model 
we show that there is similar behaviour of transfer and cascade for its energy E a and 
enstrophy Q a . We follow the techniques in [16] to arrive at this conclusion. Here we 
summarize our main results. The complete details can be found in the Appendix. 

For m ^> kf, let e ^ be the net amount of energy E a (see (3)) transferred per unit 
time into the modes higher than or equal to m, and rf^ be the net amount of enstrophy 
Q a (see (5)) transferred per unit time into the modes higher than or equal to m. From 
the scale-by-scale evolution equation of energy E a and enstrophy Q a we get the following 
relationship between and 77° 

( 18 ) ( £ m) < 27TT ZL ~2 2V 

where (•) denotes an ensemble average with respect to infinite time average measure (see 
Appendix). This result suggests that for large m, the average net of transfer of energy 
to high modes is significantly smaller than the corresponding transfer of enstrophy. This 
yields the characteristic direct enstrophy cascade. 

The inverse energy cascade is expected to take place in the range below the kf. Within 
this range, i.e. m <C kf, we obtain 

(19) (-V a J<m 2 (l + a 2 m 2 )(-e a m ), 

that is, the (inverse) average net transfer of energy to lower modes is much stronger 
than the corresponding enstrophy transfer which yields the characteristic inverse energy 
cascade. This establishes the expected 2d NS-a cascades for energy E a and enstrophy 

2.2. NS-a model effects on the scaling of the energy spectrum in the enstrophy 
cascade regime. We next derive, using semi-rigorous arguments [13, 14], the expected 
scaling for the 2d NS-a energy spectrum. We start by splitting the flow into the three 
wavenumber ranges [l,k),[k,2k),[2k,oo). We assume kf < k, where kf is the forcing 
wavenumber, since we are interested on the effects of the NS-a model in the enstrophy 
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cascade regime. We decompose the u, v and vorticity q corresponding to the three 
wavenumber range, and for simplicity, we denote as follows: 

U = U k + U k,2k + U2k 
V = + V K2 k + ^ 

Q = Qk + Qk,2k + q 2k , 

where, 4>f — 4>i,h 4>i = 0«,oo- We recall (•, •) and | ■ | denote the L 2 -inner product and 
L 2 -norm respectively. The enstrophy balance equation for the NS-a model for an eddy 



of size ~ k 1 is given by 

1 d 

2 dt 



(20) ^ 37(^,2^ qk,2k) + ^(-Ag fci2fc , qk,2k) = Z k - Z 2k , 



where, 



~Z 2k 

b(u, v, w) 



-K u k> ?* » <?fc,2fc + q 2k ) + b(u k ,2k + «2fe> <ik,2k + q 2k , qt) 

-K u 2k> ?2fc, qk + qk,2k) + b(u% + u k , 2k , q k + qk,2k, q 2k ) 
I'm • Vf , u>) . 



may be interpreted as the net amount of enstrophy per unit time that is transferred 
into wavenumbers larger than or equal to k. Similarly, Z 2k represents the net amount of 
enstrophy per unit time that is transferred into wavenumbers larger than or equal to Ik. 
Thus, Z k — Z 2k represents the net amount of enstrophy per unit time that is transferred 
into wavenumbers in the interval [k, 2k). Taking an ensemble average (with respect to 
infinite time average measure) of (20) we get 

(21) v ((-Ag fej2 fc, q k ,2k)) = (Z k ) - (Z 2k ) ■ 

Using the definition of E a (k) in (13), we can rewrite the averaged enstrophy transfer 
equation (21) as 

r2k 

uk 5 (l + a 2 k 2 )E a (k) ~v k 4 (l + a 2 k 2 )E a (k)dk ~ (Z k ) - (Z 2k ) . 

Jk 

Thus, as long as uk 5 (l + a 2 k 2 )E a {k) <C (Zk) (that is, (Z 2k ) ~ (Zk), there is no leakage of 
enstrophy due to dissipation), the wavenumber k belongs to the inertial range. Similar to 
the other a-models, the correct averaged velocity of an eddy of length size ~ 1/k is not 
known. That is, we do not know a priori in these models the exact eddy turnover time 
of an eddy of size ~ 1/k. 

In the forward cascade inertial subrange, Kraichnan theory [22] postulates that the 
eddies of size larger than 1/k transfer their energy to eddies of size smaller than 1 / (2k) 
in the time r k it takes to travel their length ~ 1/k. That is, 

(22) n ~ -L 

kU k 
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where Uf. is the average velocity of eddies of size ~ 1/k. Since there are two velocities in 
the NS-a model, there are three physically relevant possibilities for this average velocity, 
namely 

1/2 / r 2k \ 1/2 



/ 1 r \ ' f r 2k \ 

U° k = JJv k ,2k(x)\ 2 dx) ~ I J (l + a 2 k 2 )E a (k)dk) ~ (k(l + a 2 k 2 )E a (k)) l/ \ 

U l = \j2 J u k,2k( x ) ■ Vk,2k(x)dx) ~ I J E a (k)dkj ~ (kE a (k)) 1/2 , 
^ /I /"I ^,2, \ V2 f f 2k ^M_, h \ 1/2 f kE a (k) \ 1/2 
That is, 

( 23 ) U k ~ (l + a 2 fc 2)(n-l)/2 ' fa = 0, 1, 2). 

We may therefore write the turnover time scale of an eddy of size /c _1 in (22) as 

< 24) ^~w = Wi))^ ■ ( " = ' 1 ' 2 '- 

The enstrophy dissipation rate t] a which is a constant equal to the flux of enstrophy 
from wavenumber k to 2k is given by 



1 f 2h 

(25) r/ Q ~ — / (1 + a 2 k 2 )k 2 E a (k)dk 

T k Jk 



'k Jk 

and hence 



2 (E a (k)f 2 
(1 + a 2 k 2 f n -V/ 2 '' 

_ , M ^ /3 (l + a 2 P)("- 3 )/ 3 



k 3 

Thus, the energy spectrum of the smoothed velocity w is given by 

E (k) { ^ /3fc ~ 3 ' when ka < 

(26) *W - ~ when to » , 

I Q,2(6-n)/3 

Therefore, depending on the average velocity of an eddy of size k~ l for the NS-a model, we 
obtain three possible scalings of the energy spectrum, ^-( 21 - 2n )/ 3 ; ( n = o, 1, 2) all of which 
decay faster than the Kraichnan k~ 3 power law, in the subrange ka 3> 1. We note that in 
the work of [26] the k~ 17 ^ 3 power law was obtained using dimensional analysis, consistent 
with n = 2 in our notation, using a time scale which depends solely on the smoothed 
velocity field. The actual power law was not determined at the time due to insufficient 
dynamic range of their simulations. Our simulations here show that the scaling is k~ 7 , 
which corresponds to n = in our notation. 
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3. Results 

3.1. Details of the numerical simulation. The Navier-Stokes equations and the Navier- 
Stokes-a model equations with stochastic forcing were solved numerically in a periodic 
domain of length L = 1 on each side. The wavenumbers k are thus integer multiples 
of 27r. A pseudospectral code was used with fourth-order Runge-Kutta time-integration. 
Simulations were carried out with resolutions ranging from 256 2 to a maximum resolution 
of 4096 2 on the Advanced Scientific Computing QSC machine at the Los Alamos National 
Laboratory. To maximize the enstrophy inertial subrange, the forcing is applied in the 
wavenumber shells 2 < k < 4. We also add a hypoviscous term (fiA~ 1 v) which provides a 
sink in the low wavenumbers. Let 7i — 1 — a 2 A, then v = Tiu. We simulate the following 
equation: 

(27) d t u - n~\u x (V x v)) = -Vtt + fx(A~ l u) + vAu + / 

where 7r denotes the modified pressure and /i > is a hypoviscosity coefficient. The 
equation in (27) is equivalent to 

(28) d t v - u x (V x v) = -Vp + niA^v) + uAv + /, 

where / = TC(f) (compare with equation (2)). We demonstrate in the next section, that 
adding a hypoviscous term does not have a significant effect on our observed numerical 
results. We use no extra dissipation term in the small scales besides the regular dissipation 
in the Navier-Stokes or the Navier-Stokes-a. Some relevant parameters and results of our 
simulation are presented in Table 1. The 256 2 simulations were performed to ascertain 
whether it was feasible to perform the computations in reasonable time without a low- 
wavenumber sink term, the hypoviscosity. The simulations of 1024 2 and higher show the 
convergence of the scaling exponent of interest, namely the value of 7. 

The condition for a resolved simulation is judged by the degree to which the enstrophy 
dissipation scale is resolved since the enstrophy cascade is expected to govern the dynamics 
in the range k ^> kf. The enstrophy dissipation rate for NSE and NS-a are: 

(29) rj = rj a = u(\Vq\ 2 ) 

where q = V x v, and v is respectively, either the NS velocity or the unfiltered velocity for 
the NS-a. The enstrophy dissipation scale is then computed, following Kraichnan [22, 23], 
by dimensional analysis in a manner analogous to the way in which the Kolmogorov energy 
dissipation scale is calculated in 3d turbulence: 

/ 3 \ V6 

(30) i, M =(^j ■ 

A resolved simulation has fc max /^ > 1. In all our simulations, k max l v > 2, indicating 
well- resolved flow. Note, however, in table 1 that, keeping all other parameters fixed, 
increasing a decreases k max l v , indicating that as if the NS-a flow is less resolved, from 
the point of view of the enstrophy cascade, than the NSE for a given grid and viscosity. 
However, this observation is misleading since the computations for the NS-a were done 
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Table 1. Parameters of the simulations: number of grid points per side N, a in units of 
Ax, viscocity coefficient hypoviscosity coefficient /x, maximum wavenumber k max = 
V2 ... /^\ 1/6 

— N , enstrophy dissipation scale l v = — , where r\ a is the enstrophy dissipation 
rate. 
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with a rescaled forcing term with an amplitude growing indefinitely as a-^oo. Therefore, 
by increasing a, the NS-a is forced more vigorously, and damped strongly. We discuss 
this further in section (3.4). 

3.2. Dependence of scaling behaviour on hypoviscosity. The hypoviscosity term 
in (28) provides a sink in the low wavenumbers for the energy. This allows the flow to 
reach statistical equilibrium in more reasonable computational time. However, since it is 
an ad hoc addition to the NSE or NS-a model, we need to ascertain whether it affects 
the behaviour in the range of scales of interest, namely the range k > 1/a. 

In figure 1 we compare the energy spectra for 256 2 run for DNS (a = 0) and 2d 
NS-a with a = 2.04 (in units of Ax). The solid lines represent the DNS runs with 
varying hypoviscosity coefficient. In the inset, the compensated plots for the DNS show 
very little difference. The dotted lines correspond to the 2d NS-a runs with varying 
hypoviscosity coefficient. Again, as seen in the compensated plots in the inset, there is 
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Figure 1. Energy spectra for a 256 2 simulation with fixed viscosity and varying hypo- 
viscosity coefficient fi. The wavenumber k is in multiples of 27r. The solid lines are the 
DNS (a = 0) calculations of E(k). The dotted lines are the NS-a model calculations of 
E u {k) for small a. The behaviour of the spectra is largely independent of the magnitude 
of the hypoviscosity in the enstrophy cascade subrange (6 < k < 15). The inset shows 
the spectra compensated by fc 4 5 . The resolution of this simulation is far to small to 
observe the expected scaling exponent. 



very little dependence on the hypoviscosity. These numerical results show that for a small 
to moderate hypoviscosity coefficient, the spectral slope in the enstrophy inertial range 
is not significantly affected by the addition of the hypoviscous term. Furthermore, as we 
see in Table 1 for the 256 2 simulation, the enstrophy dissipation length scale l v is not 
affected by varying hypoviscosity coefficient. 

From the above empirical observation, we conclude that, for the range of scales of 
interest, we can safely use a hypoviscous term and save significant computational time in 
our higher resolution runs. 

3.3. Varying alpha; convergence to NS-oo. In this section we show the effect of 
varying the parameter a. For a given grid, we fix the viscosity coefficient and vary a in 
order to measure the scaling exponent 7 of the energy spectrum in the enstrophy inertial 
subrange. 

In figures 2, 3 and 4 we show the NSE spectrum E(k) and the spectrum E u (k) from 
1024 2 , 2048 2 and 4096 2 simulations, each with fixed viscosity and varying a. As expected, 
the scaling ranges increase as the number of grid-points increase. In each figure the solid 
black line is the DNS spectrum E(k), and approaches a scaling close to k~ 3 as N increases. 
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k 



Figure 2. Energy spectra for 1024 2 simulation. The black curve is the DNS (a = 0) 
which shows close to k~ 3 scaling in the cnstrophy cascade range 6 < k < 20. The solid 
red curve is the E u (k) spectrum as a — > oo which scales close to k~ 7 in the enstrophy 
cascade range 6 < k < 25. The energy spectra for intermediate values of a tend to 
the a — -> oo limit as a increases. The inset shows the DNS energy spectrum (black) 
compensated by k 3 ' 7 and the a — > oo energy spectrum (red) compensated by k 7A 

In both Figs. 2 and 3 for < a < 15 (in units of Ax), we see the spectrum E u (k) peels 
away from the NSE spectrum (a = 0) near ka — 1 but displays no clear scaling behaviour 
for ka > 1 until a > 100. To discern a clear power-law of the NS-a model spectrum, we 
consider the data from simulation of the NS-oo equations (17). This allows us to see the 
scaling of the NS-a model energy spectrum without contamination by finite-a or DNS for 
the NSE (k~ 3 ) effects. At our maximum resolution of 4096 2 in Fig. 4, as a— >oo (solid red 
line), there is a clear convergence of the scaling to k~ 7 . Table 2 summarizes our findings 
for the scaling of the E u (k) spectrum as a— >oo for each of our simulations. According to 

Table 2. Convergence of a— >oo scaling exponent 7 of E u (k) as the resolution is increased. 



N 256 1024 2048 4096 
7 8.0 7.4 7.1 7.0 



(26), the scaling k 7 corresponds to an enstrophy turnover time scale determined by the 
velocity v. We conclude that the dynamics of the smoothed velocity field, which is the 
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Figure 3. Energy spectra for 2048 2 simulation. The wavenumber is in multiples of 2ir. 
The black curve is the energy spectrum of the DNS which shows close to k~ 3 scaling 
in the enstrophy cascade range 6 < k < 35. The solid red curve is the E u (k) spectrum 
as a — > oo which scales approximately as k~ 7A in the wavenumber region 6 < k < 25. 
The inset shows the DNS energy spectrum (black) compensated by fc 3 5 and the a — > oo 
energy spectrum (red) compensated by fc 71 

puported model for turbulence, is nevertheless still governed by the unfiltered velocity 
field. 

In this analysis and conclusion, we have assumed that the scaling in the asymptotic 
limit a^oo would be identical to the scaling for ka > 1 for small a. In Fig. 5 we 
demonstrate that this assumption is reasonable by showing that the scaling of the energy 
spectrum for small a = 6.5 is approaching close to k~ 7 scaling for a small range of ka > 1. 
We are therefore reasonably convinced that going for the asymptotic limit is indeed giving 
us the correct scaling for finite (small) a; the a^oo limit merely maximizes the range 
over which a clear scaling exponent can be measured at a given resolution. 

3.4. NS-q model effects on the dissipation length scales of the flow. In Table 1, 
we observe that, keeping the viscosity, hypoviscosity and the resolution fixed, increasing 
a tends to decrease the enstrophy dissipation length-scale l v (see (30)). We also learned 
that the rate of dissipation of the energy E u is dictated by the eddy turnovertime of the 
unfiltered velocity field. In this section, we further explore the effects of the a parameter 
on the smallest scales of the flow. We recall that the NS-a computations were done with 
a rescaled forcing term, a fact which is important in the discussion to follow. 



14 




10° 10 1 10 2 10 3 

k 

Figure 4. Energy spectra for 4096 2 simulation. The black curve is the spectrum for the 
DNS, the red curve is the spectrum for a — > oo. The black curve in the inset corresponds 
to the NSE energy spectrum compensated by fc 3 ln(fc/ + fc) 1 / 3 following [23]. The red 
curve in the inset is the energy spectrum E u (k) for NS-oo compensated by k 7 . The 
region 6 < k < 40 is flat indicating the nominal range over which the k~ 7 scaling holds. 

Table 3. Dissipation length scales when varying a. N = resolution, a, v = visocisty 
coefficient, fj. - hypoviscosity coefficient, l v - enstrophy (O a ) dissipation length scale, l Vu - 
smoothed enstrophy dissipation length scale 



N 


a (in units of Ax) v 






N 


a a v 


/' 




1024 


le~ 4 


15 


.004133 


1024 


le" 4 


15 


.004133 




3.25 




.004099 




3.25 




.004490 




15 




.003973 




15 




.005659 




100 




.002165 




100 




.006827 




30 




.000488 




3C 




.006858 



In the left panel of table 3, we present the enstrophy dissipation length scale l v (see 
(30)) for the 1024 2 simulation as a increases. This length scale, as we already know, gets 
smaller with increasing a value. A visual of this effect is given in Fig. 6, where we plot 
the isosurfaces of vorticity for increasing values of a. Observe how the vorticity values 
grow while the vorticity structures become much finer as a increases. 
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Figure 5. Compensated energy spectra for 2048 2 simulation for a = 6.5 (k a = 39.75; 
vertical dashed line). The energy spectrum is compensated by k 7 ,k 19 ^ 3 , and , k 17 ^ 3 
respectively. The region 39 < k < 70 in the first subplot follows a flat regime which 
indicates the nominal range over which the k~ 7 scaling holds. 

By contrast, consider the right panel of table 3 where we present a "dissipation" length 
scale 

(3D , 

(32) where i] u = v 

(33) and u = V x u. 

corresponding to a naively calculated length-scale for the smooth enstrophy Q u = |Vxit| 2 . 
This length-scale grows as a increases. Corresponding to this we present in Fig. 7 the 
isosurfaces of V x u for increasing values of a. Note that in this case, the vorticity values 
diminish while the vorticity structures become increasingly smooth and diffuse. 

Thus, on the one hand the smooth velocity field and its vorticity are consistent with 
reduced resolution requirements for the NS-a model; on the other, the behaviour of the 
conserved quantity for the (inviscid and unforced) NS-a model indicate a requirement 
for increased number of grid points, counter to what one would like to see in a sub-grid 
model. Note, however, that in order to avoid trivial dynamics as a increases, we have 
scaled the forcing term in the simulations of the NS-a. That is, the computed Eq. (28) 
tends to the Eq. (17) as a — > oo. This was done so that we could conveniently study 
the case of large a, thus extending the scaling range of interest. It could well be that for 




Figure 6. Isosurfaces of vorticity V x v for the 1024 2 simulation, a = 0, 3.25, 15, 100, oo, 
reading each row of figures from left to right. The vorticity field exhibits increasingly 
fine structures as a is increased. 




Figure 7. Isosurfaces of vorticity V x u for the 1024 2 simulation. a = 
0,3,25,15,100,oo, reading each row of figures from left to right. The structures be- 
come smoother with increasing a. 



IT 



a detailed study of small-a with the forcing unsealed, the desired computational gains 
expected of a sub-grid model could be observed. However, our main goal in this study, the 
scaling exponent 7, could not be clearly obtained at achievable resolutions without going 
for the large-a limit. Thus, a straightforward comparison, for the purpose of checking 
the NS-a as a sub-grid model, is not clear. This implies that the implementation of the 
NS-a model as a turbulence model needs to be performed with some care. 
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Appendix A. Average transfer and cascade of the energy E a and 

ENSTROPHY fi Q IN THE TWO-DIMENSIONAL NS-« MODEL 

Here we will show analytically, the transfer and cascade of the energy in the two- 
dimensional NS-a model. Following the exposition as in [16], we will show that the 
conserved (in the absence of viscosity and forcing) energy E a (3) and enstrophy Q a (5) 
in the two-dimensional NS-a model have similar transfer and cascade behaviour as in the 
two-dimensional turbulence. We recall that we denote by <pf = <fr lt i and <pf = 0; jOO 

(1) Average energy transfer 

We start by decomposing the velocity fields into two components: the large-scale 
component w< and v< containing eddies of size larger than 1/m and the small- 
scale component w> and u> containing eddies of size smaller than or equal to the 
lengthscale 1/m. That is, 



u — u m + u m 
v = v < + v > . 



For the forcing /, we assume it to be time independent and contains only finite 
number of modes, namely, 

/ = fm,m, where, 1 < m < m < 00. 

Here we will look at the two cases: m e (1, m) and m £ (m, 00). We recall that 
(•, •) and I • |, denote the L 2 -inner product and L 2 -norm, respectively. From the 
above decomposition, we define the kinetic energy contained in the large and small 



eddies (respectively) as 



E a = \ («m> V m) alld E * = \ ( u m» V m) ■ 

We are interested on how the energy i£ Q is transferred between large and small 
scales. Consider the case where m > m (i.e. beyond the injection of energy). 
Following the notation as in [14], let us denote by b(u, v, w) = {—u x (V x v), w). 
We write (2) in its functional differential form (see,e.g. [16, 28]), and apply the 
large and small scale decomposition to get the corresponding energy equations for 
the large and small scales 



te U m) + » Om; ~ Au m) + K U m + U m, V m + V m> U m) = if, U m, 



d 

dt 



Define 

$< : = -b(u>, v>, u<) + b(u< , v<, w>) 
$> : = b(u> , v>, u<) - b(u<, v<, «>) 

Notice that $ > + $ < = 0. We can rewrite the energy equations in the following 
form 

j t E< + v (v<, -Au<) = $< + (/, u<) 
j t E> + u(v>,-Au>) = <t>>. 

We can interpret the individual terms above as follows: (/, w< ) - represents the 
energy flow injected into the large scales, by the forcing term. $ < - represents the 
net amount of energy per unit time that is transferred from small to large length 
scales. $ > - represents the net amount of energy per unit time that is transferred 
from large to small length scales. — b(u^, v<, it>) - represents the energy flow 
induced in the high modes by inertial forces associated with lower modes. 

Let T be the set of all vector trigonometric polynomials with periodic domain 
Q. We then set V = {</> G T : V • <j> = and f n <j>{x) dx = 0}. We set V to be the 
closure of V in the Sobolev space H 1 . The solution operator u{t) = S a (t)u of (2) 
defines a dynamical system on the phase space V. A generalized notion of infinite 
time averaging, as it was defined in [16], induces a probability invariant measure 
with respect to S a (t). We denote this measure by dP, and the ensemble average 
with respect to this measure, will be denoted by 

(0(-)> = / <j>(w)dP(w). 
Jv 
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Thus, beyond the range where energy is injected, the average transfer of energy is 
into the higher modes, that is 

) <*>)=!/<«>, -AO >0. 

Now consider modes that are below the injection of energy For that purpose, 
assume m > 1 and consider m such that 1 < m < m. The NS-a equation (2) can 
be decomposed as follows: 

J t («m> U m) + v ("m. ~Au>) + + U> , W< + V>, U>) = (/, M> ) . 

The associated energy equation for the lower modes reads: 

|^ + r/fc-A«<) = $< 

We then take the ensemble average to get 

($<> = v (v<, -Au< ) > 0. 

That is, for wavenumber regime below the injection of energy the average transfer 
of energy is from high wavenumbers to lower wavenumbers. 

(2) Average enstrophy transfer 

Next we consider the details of the transfer of enstrophy. Again, we assume / 
is time independent and contains only a finite number of modes. Here, we do a 
slightly different exposition as in [16]. We take the curl of the momentum equation 
and using the incompressibility condition and the vector identity, 

V x (a x b) = a(V • b) - 6(V ■ a) - (a ■ V)6 + (b ■ V)a, 

we get the vorticity formulation for the 2d NS-a equation: 

d t q — uAq + u ■ Vq = V x / 

where, q = V x v and we recall that v = u — a 2 Au. We split the flow into two 
components u = w< + w> and q = q< + g>. 

The amount of enstrophy contained in the large and small eddies are given, 
respectively, by 

^ = |Vx^<| 2 
^ = |Vx^>| 2 . 

We would like to show how the enstrophy is transferred between the large and 
small scales. First consider the case m > m (i.e. beyond the injection of energy). 



Following the notation as in the NSE [16], we denote by b(u,v,w) = (u ■ Vv,w). 
We write the evolution of the enstrophy governing the large and small scales: 

J t (?m> Qm) + V {Qmi ~ A Qm) + %m + «m- Qm + 9m. 9m) = (V X /, <?<) 
^ (?m> ?m) + V (<lm> ~ A Qm) + h ( U m + U>, q< + q>) = 0. 

where, we denote by 

^ • b{ U m ~l~ U m) 9m ~l~ 9m; Qm) ^{ U mi Qmi Qm) ^{ U mi Qmi Qm) 

^ ■ b{ U m U mi Qm Qmi Qm) ^i U mi Qmi Qm) ^( U mi Qmi Qm) 

\l/ < - is the net amount of enstrophy transferred into low modes and \E f> - is the 
net amount of enstrophy transferred into high modes. Note that for almost every 
t, \I /< + ^> > = 0. From this we can rewrite the enstrophy equations above in the 
following form 

^kml 2 + HVg<| 2 = v[/< + (Vx/,g<) 



2dt 

2dt \Qm\ ^^ivy m l 

Taking the ensemble average (with respect to the infinite time averaging measure 
dP) of the equation for enstrophy associated with the low modes, we get 

<vf>>=^|Vg>| 2 )>0. 

This implies that, beyond the injection of energy, the average transfer of enstrophy 
is from low modes into higher modes. Next we consider the modes below the 
injection of energy. Assume 1 < m < m, we get 

j t (Qmi Qm) + v (<?<, -Ag<) + b(u< + u>, q< + q>, q<) = 

j t (Qmi Qm) + v (q m , -Ag>) + b(u< + u> , q< + q>, q>) = (Vx/, q>) . 

We take the ensemble average of enstrophy equation associated with the low modes 
to get 

<*<> = -<*>> = HVg<| 2 >0, 

that is, the average net transfer of enstrophy is from high modes to lower modes 
for wavenumber regime below the injection of energy. 

(3) Direct enstrophy cascade and inverse energy cascade 

In (i) and (ii) we have seen that in the range above the injection of energy both 
the energy and enstrophy is transferred from low to higher wavenumbers. On the 
wavenumber s regime below the injection of energy, both the energy and enstrophy 
are transferred from high to lower wavenumbers. Here we will show that, in the 



21 



\Qm',m"\ +^|Vg 



|2 i ,.|V7„ |2 

lm'm" 



wavenumbers regime above the injection of energy, there is a much stronger transfer 
of enstrophy leading to what is called direct enstrophy cascade . On the other hand, 
in the wavenumbers regime below the injection of energy, the energy transfer is 
much stronger than the enstrophy transfer leading to what is known as the inverse 
cascade of energy. 

We split the flow into three parts assuming the same form restriction on the 
forcing as above. Consider m" > m' > m. Let u = u m , + w m ',m" + u m „ and 
Q = Qm' + Qm',m" + Q m " ■ The vorticity field associated with the wavenumbers 
between vnl and vnl' is 

^?m',m" — vAq m i tm ii + B{u m , + u m > :7n » + u m „,q m , + q m > , m " + Q m ") = 

where we denote by B(v , w) := (v ■ V)w. The evolution equation for the enstrophy 
associated with the modes between m' and m" is given by 

ld_ 

2dt' 

= ~K U m> + u m',m" + U m ',, q m , + q m \m" + Q m " , 1m' ,m") ; 

where, b(u, v,w) := ((«• V)f , w). One can compute by using the bilinear properties 
of b (see, e.g. [11, 28, 16]) that, 

-b(u m , + u m , tTn n + u>„, q m , + q m ,^ m „ + q>„, q k ) = 7] m , - r] m „ 

where, 

r] m , - is the net amount of enstrophy transferred per unit time into the modes 
higher than or equal to m' and r] m „ - is the net amount of enstrophy transferred 
per unit time into the modes higher than or equal to m", and given by 

Vm' — ~K U m' > Qm' ' Qm' ,m" + Qm") + K U m' ,m" + U m „, q m ' ,m" + Q m » , Q m >) 
~Vm" = ~K U m" ) Qm" 1 Qm' + Qm' ,m") + K U m' + U m' ,m" , Q m < + Qm' ,m" j Q m >) ■ 

In explicit form we have 

Take the ensemble average of (50), we get 

(Vm") = (rjm') -^{\^Qm',m"\ 2 ) ■ 

That is, the average net transfer of enstrophy into the modes higher than m" 
is equal to the average net transfer of enstrophy into the modes higher than or 
equal to vnl minus the enstrophy lost due to viscous dissipation within the range 
[m',m"]. In [m',m"], it is assumed that the viscous dissipation is negligible and 
the enstrophy is simply transferred to smaller and smaller eddies. This occurs 
whenever 

(Vm")>(Vm'}^>^(\^Qm',m"\ 2 ) 



The range of wavenumbers k > m up to where (51) holds is called the enstrophy 
inertial subrange. Now within this range there is still an average net transfer of 
energy to higher modes. Denote by (e^) : = ($ > ) and (77^) := (^ > ). From (38) 
and (47) 

1 , « v < (\qm,oo\ 2 ) (I Vg mj00 | 2 ) _ 1 (C) 

\ £ ml — 1 , „.<> 2 — 



v m 1 + a 2 m 2 m 2 (l + a 2 m 2 ) v m 2 (l + a 2 m 2 ) 
Hence, 



m 2 (l + a 2 m 2 ) 



This result suggests that for large m, in particular, kf < m < kj, where kj is the 
dissipation wavenumber, the average net of transfer of energy to high modes is 
significantly smaller than the corresponding transfer of enstrophy. This yields the 
characteristic direct enstrophy cascade in this range. 

The inverse energy cascade takes place in the range below the injection of energy. 
Consider 1 < m! < m" < m. We follow the same steps as above. We decompose 
the flow into three components u = + u m ^ m n + and proceeding as before 
we obtain 

(e£„) = (e™,) + v ((vm'm'', -Au m > >m »)) 
As we have seen in (i) we have in this case the inverse transfer of energy 

(e«„) < and (e«,> < 

Therefore, as long as 

( £ m") ^ (Ci') < V (Vm>,m", -A'U m / jm ») 

we have inverse cascade. Now within the range corresponding to energy cascade 
(below injection of the force /) both the energy and enstrophy are transferred to 
lower modes 

(-O = v {(Vv hm , V«i, m )> = v (I Vn 1/m | 2 + o?\ Au ljm | 2 ) > 
(-C) = ^<|Vg lim | 2 >>0 
Since U\^ m contains only modes smaller than m, we therefore have 
I Vgi, m | 2 < m 2 (l + a 2 m 2 ) (Vui, m , Vui, m ) 

which implies 

(-C><™ 2 (i + « 2 ™ 2 M-C> 

that is, in the wavenumber regime below the injection of energy, the (inverse) aver- 
age net transfer of energy to lower modes is much stronger than the corresponding 
enstrophy transfer which yields the characteristic inverse energy cascade. 
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